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ABSTRACT 


Multirate digital signal processing techniques have been developed in the re¬ 
cent years for a wide range of applications, such as speech and image compression, 
digital audio, statistical and adaptive signal processing, numerical solution of differ¬ 
ential equations and many other fields. 

The purpose of this thesis is to extend optimal filtering techniques to random 
signals sampled at different rates. In particular, two major problems are considered: 
(1) optimal filtering of two sets of observations at different sampling rates as a multi¬ 
rate Wiener filter, and (2) linear prediction on successive samples of a random process. 
In the first problem it is shown that the standard Wiener filter can be extended to the 
multirate case, while preserving its optimality. In the second problem it is shown that 
multichannel linear prediction on successive samples of a process, yields orthogonal 
uncorrelated innovations. 
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DISCLAIMER 


The computer programs in the Appendix are supplied on an “as is” basis, with 
no warrantees of any kind. The author bears no responsibility for any consequences 
of using these programs. 
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I. 


INTRODUCTION 


In many practical applications of digital signal processing we have to process 
signals sampled at different rates. For example, in telecommunication systems that 
transmit and receive various types of signals (e.g. teletype, facsimile, speech, video, 
etc.), the different sampling rates are related to the bandwidths of the corresponding 
signals. The process of converting a signal from a given sampling rate to a different one 
is called sampling rate conversion. In turn, systems that process signals at different 
rates are called multirate digital signal processing systems. 

The most practical method of sampling rate conversion of a digital signal is 
to perform it entirely in the digital domain, as a combination of upsampling and 
downsampling by an integer factor [Ref. 3]. In particular, the process of reducing the 
sampling rate by an integer factor D (downsampling by D) is called decimation, while 
the process of increasing the sampling rate by an integer factor U is called upsampling 
(it is also called expansion or interpolation). 

A. PREVIOUS RELATED RESEARCH 

During the past decade a very important area of research has been in Mul¬ 
tiresolution and Multirate Digital Signal Processing. Two very successful approaches 
to this problem have been developed: one based on wavelet decomposition [Ref. 1], 
and another based on the multiresolution Kalman filtering [Ref. 4] and [Ref. 5]. In 
the following sections, we provide an outline of both methods and their applications 
in multirate signal processing. 

1. Wavelet Decomposition 

A wavelet is a “small wave,” a wave that has its energy concentrated in time 
to give a tool for the analysis of transient, non stationary, or time-varying phenom¬ 
ena [Ref. 1]. The property stated above, allows the wavelet to perform analysis in 
both frequency and time with a flexible mathematical foundation. For that purpose, 
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the wavelet decomposition uses two functions; the scaling function ip(t) that is respon¬ 
sible to track the signal in time, and the mother wavelet ip(t), related to (p(t), that is 
the function that performs analysis of a signal using multiple resolutions. Based on 
the two functions defined above, a large number of signals can be decomposed as an 
expansion also called the Discrete Wavelet Transform (DWT): 

OO OO OO 

/(<)= £ <*•*>(«-*)+ £ (i.i) 

k=—oc k=—ooj=0 

The functionality of Equation 1.1 can be described with the use of vector 
spaces. We call S the vector space of signals where any f(t) € S can be expressed as 

f(t) = 52a k -<p k (t). (1.2) 

k 

Then, the set of functions (p k (t ) is called an expansion set for S. If the representation 
is unique, the set is a basis. If one starts with the basis set and defines the space 
S as the set of all functions /(f) that can be expressed by Equation 1.2, this set is 
called the span of the basis set. Then, we define a set of scaling functions in terms of 
integer translation of the basic scaling function by 

<Pk(t) = <p(t - k) (1.3) 

where k is an integer (k 6 Z). One can generally increase the size of the subspace 
spanned by changing the time scale of the basis functions. However, the important 
features of a signal can be described or parameterized more efficiently by defining a 
set of functions 

Kk{t) = 2^{2H - k) (1.4) 

that span the “difference” subspaces or orthogonal complement spaces, defined as the 
subspaces between the spaces spanned by the various scales of the scaling function. 

The selection of the mother wavelet is not unique, however tp(t) is usually 
chosen such that the set of functions in all resolutions comprises an orthonormal set , 
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where for a complex basis set: 

(^fc(*)»^(*)> = / ^k(t) • t*(t)dt = 0 (1.5) 

and 

= i- (i-6) 

The need for wavelets came from the inherent weakness of the Fourier-based methods 
to analyze signals with edges, time-varying signals, or broadband signals such as 
transients. On the other hand, the DWT uses a 2-D expansion set (the indexes j for 
frequency and k for time) that provides a very close to optimal analysis of a signal, 
localized in both time and frequency. 

A particularly interesting class of problems of multirate signal processing is 
multiresolution analysis, where a signal is decomposed into a lower and an upper 
frequency component in a recursive fashion [Ref. 6: pp.254-259]. The filter banks 
that implement this kind of structure are called tree structured filter banks and they 
are very useful in the implementation of wavelet analysis. In particular, one can 
express xf(t) and ip{t) in terms of filtering functions as follows: 

(pit) = JZ h 0 (n) • V2 • <p(2t - n) (1.7) 

n 

where n is an integer (n € Z), and the coefficients ho{n) are a sequence of real or 
complex numbers called the scahng function coefficients (or the scaling filter) and y/2 
maintains the norm of the scaling function. Equation 1.7 is called the multiresolution 
analysis (MRA) equation, or the dilation equation. Also 

ip(t) = hi(n) ■ y/2 • <p(2t - n) (1.8) 

n 

for some set of coefficients hi(n). From the requirement that ipj,k{t) span the “dif¬ 
ference” subspaces, and the orthogonality of integer translates of the wavelet, the 
coefficients h\{n) are required to be related to the scaling function coefficients by 

h 1 (n) = (-l) n -h 0 {l-n). (1.9) 
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We can then perform an analysis from fine scale to coarse scale (fine to coarse resolu¬ 
tion) and obtain the coefficients c*, and dj^ that appear in Equation 1.1 by iterating 
the tree-structured sections of Figure 1. 



Figure 1. Two-Stage Two-Band Analysis Tree (from [Ref. 1]) 

Despite the advantages that were previously stated, the wavelet theory can 
only perform sub-optimal analysis of a signal, whether it is stationary or not. The 
reason for this is that it does not take into account the signal properties or, in the case 
of a random signal, the signal statistics. In addition, wavelet theory is specialized in 
the analysis and synthesis of individual signals, whether they are 1-D or 2-D. However, 
there is a new class of applications that deals with the synthesis of different signals 
that descent from the same origin and collected individually (i.e. different sensors), 
called information fusion. This is the kind of problems that the following application 
attempts to solve. 
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2. Multiresolution Kalman Filter 

Multiresolution Multirate techniques have been developed during the past ten 
years, based on the Kalman Filtering approach. The result is a recursive algorithm 
both in time as well as in resolution, which optimally combines estimates and obser¬ 
vations at different levels of resolution and sampling rate. The general theory has 
been presented initially in [Ref. 7], and then refined in [Ref. 4] and [Ref. 5]. The 
framework presented in this fine of research is based on a general solution of a Riccati 
equation, recursive in time and in resolution. 

A different approach has been taken by [Ref. 8], where multiresolution mod¬ 
eling and filtering of time series have been addressed based on a recursive Kalman 
Filtering approach. It is shown that a time series can be decomposed into a number 
of models operating at different sampling rates, having independent innovations. Pre¬ 
diction and optimal filtering based on a number of observations at different sampling 
rates have been presented. 

The general Kalman Filtering approach yields a class of estimators which are 
Infinite Impulse Response (HR), and attempt to determine the optimal prediction 
of the time series at time n + 1, based on all the observations from the initial time 
(say n = 0) to the current observation time n. We say that this approach has an 
infinite memory, and therefore it necessitates HR filters to be implemented. However 
the tendency in applied Digital Signal Processing is to use Finite Impulse Response 
(FIR) filters, and the approach of this thesis will be to develop a class of optimal 
predictors based on a finite window of data. 

B. THESIS APPROACH 

While the previously mentioned approaches have been successful, they are 
somewhat difficult to interpret in terms of conventional sampling and filtering op¬ 
erations that are the mainstream of signal processing for deterministic signals. The 
arguments presented in this thesis are based on a number of methods representing 
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decimation and linear filtering. It turns out that a very effective way of representing 
them is by using linear algebra techniques. In this way each operator is associated 
with a matrix operation, and the analysis is based on a number of concepts available 
in the linear algebra literature. 

In the first part of this research we formulate the tools for developing the op¬ 
timal filters for problems with multiple observations at different sampling rates. In 
particular, we show that a Wide Sense Stationary (WSS) stochastic process can be 
estimated based on observations at different sampling rates, all affected by measure¬ 
ment noise. The outcome is an optimal linear time-varying FIR filter with periodically 
changing coefficients. The filter is designed on the basis of a set of equations similar 
to the Wiener-Hopf equations in the standard Wiener filter setting, derived directly 
from the orthogonality principle. The resulting formulation is an information fusion 
application, developed for the case of two observation sequences at different rates, 
but extendable to more than two sets of observations. 

In the second part of the research we use the methods of linear prediction and 
its multichannel extension in order to perform linear prediction of a random process 
in multiple resolutions. In particular, we form two separate channels from the even 
and odd samples of the process, and perform linear prediction on the multichan¬ 
nel process consisting of these two sets of samples. The resulting formulation is a 
multiple resolution application, extending to capabilities similar to those of wavelet 
decomposition. 

C. THESIS OUTLINE 

This thesis report consists of five chapters including the introduction. Chap¬ 
ter II provides an introduction to notation and algebraic manipulations that are used 
extensively in the work. Chapter III develops the multirate Wiener filter and per¬ 
forms performance comparison with the single rate Wiener filter. Chapter III also 
provides variations and generalizations of the multirate Wiener filter and performs 
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performance analysis simulations between the two filters. Chapter IV deals with the 
multichannel and multirate linear prediction. Finally, Chapter V summarizes the 
work and suggests future directions of research. 
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II. MATHEMATICAL PRELIMINARIES 


Wiener filtering [Ref. 2: pp.347-354], named after Norbert Wiener (1894- 
1964), is a general type of optimal linear filtering. It involves the estimation of a 
wide-sense stationary (WSS) random process s[n], called the “desired” process, that 
cannot be observed directly, from the observation of a jointly stationary random 
process x[n}. The desired random process may represent a signal which is subject 
to various forms of distortion and interference. In the simplest form of the problem, 
the goal of the signal processing is to estimate the sequence s[n] with a sequence 
s[n] from the observed sequence x[n] via a linear finite impulse response (FIR) filter. 
Specifically, it is required to use a linear combination of the present value x[n] and 
the last P — 1 values to produce an estimate that minimizes the mean-square error 
(MSE). Wiener filtering can also be performed using an infinite impulse response 
(HR) filter, however this case will not be considered here. 

In order to minimize the MSE, it is necessary and sufficient to satisfy the 
orthogonality principle. Specifically, this principle states that if we define e[n] = 
s[n] — s[n] the error in the estimation, the resulting error should be orthogonal to the 
observations, or 

E {x[n — i\ • e*[n]} = 0, where i = 0,1,... ,P — 1. (H.l) 

The filter impulse response that satisfies the orthogonality principle is denoted by the 
vector h and is obtained as the solution of the matrix equation 

Rxi • h = t sx (II.2) 

where R xx is a Toeplitz matrix formed from the autocorrelation function (ACF) of the 
sequence x[n] and r sx is a vector formed from the cross-correlation between s[n] and 
x[n]. Equation II.2 is known as the Wiener-Hopf equation. The resulting estimate 
can be written as 

s[n] = x* T • h (II.3) 
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while the minimum MSE can be expressed as 

o* = fl„[0] - il ■ h- = B„[0] - h' T • f„ (11*4) 

where the “(tilde) operator denotes the reversed of the vector (see Section A below). 

The first problem addressed in this thesis, extends the FIR Wiener filter to 
the case of multirate observation sequences. However, the guiding principles are the 
same as those discussed above. 


A. BASIC DEFINITIONS 
1. Observation Vectors 


In the analysis of random signals it is convenient to represent finite length 
sequences as vectors. We can represent a signal s[k] on the interval n — N +1 < k < n 
with the vector 


s n 


s[n — N + 1] s[n — N + 2] 



(II-5) 


and identify it as a random observation vector for s[k}. Since the values of the signal 
are random variables and are possibly complex-valued, the vector elements may also 
be complex-valued random variables. 

It is evident from the definition that the length of an observation vector may 
vary according to our will. In addition, when we have multiple observations, not 
necessarily of the same length, say from x[n] and y[n], we can define an augmented 
observation vector as a combination of x n and y n , as follows: 


X n 


Ym 


x[n — N 4- 1] x[n — N + 2] ... x[n] j 
y[m — M + 1] y[m — M + 2] ... y[m ] 


T 


(II 6) 


z 


X n 

Ym 
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2. Reversal Operation 

In the analysis that follows, it is often necessary to reverse the order of the 
points in a given vector by turning the vector upside-down (see [Ref. 2]). Given a 
vector x n we define its time reversal x n as: 



x{n — N + 1] 



x[n] 

Xn = 

x[n — N + 2] 

, then x n = 

x[n 

- N + 2} 


x[n\ 


x[n 

-N + 1] _ 


(n.7) 


Likewise, the reversal of a matrix is defined as a matrix with the elements reversed 
about both its vertical and its horizontal axis. Table I shows some elementary prop¬ 
erties of the reversal. 


- 

Quantity 

Reversal 

Matrix product 

AB 

A B 

Matrix inverse 

A" 1 

(A)" 1 

Matrix conjugate 

A* 

A* 

Matrix transpose 

A T 

A T 


Table I. Properties of Reversal (from [Ref. 2]) 


3. Correlation Matrices 

In this thesis we focus our attention on random processes. In the case that the 
processes are Gaussian in nature, the first and second moments are sufficient to fully 
characterize the process. In general, even when the processes under consideration are 
not Gaussian, the second order moments are still at the basis of the design of optimal 
linear (Wiener) filters. In particular, we define the mean of a random process as: 
m x [n] = E {z[n]}, and the autocorrelation function (ACF) of a random process as: 

Rxx[ni, n 0 ] = E {z[ni] • x*[n 0 ]}. (II.8) 
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In the particular case where the process is wide-sense stationary (WSS), the mean is 
a constant m x and the ACF is a function of the time lag between the samples: 

R xx [l] = E {x[n} ■ x*[n - l}} . (II.9) 

Similarly we define the covariance function as: C xx [l ] = R xx [l] — \m x \ 2 . When a 
random process has zero mean, the ACF and the covariance function are the same. 

Given an observation vector s n , we define the correlation matrix of a zero- 
mean, (WSS) signal s[n] as 


R SS =E{s n -s* n T } 


s[n — N 4 - 1 ] 
s[n — N + 2] 


s*[n-N + 1] s*[n — A + 2] ••• s*[n] I > • 


( 11 . 10 ) 


Due to the WSS assumption and the definition of the ACF, we can write R ss as 


Rssl o] Rss[- 1] ••• Rss[-N + 1] 

Rss[ 1 ] Rssl 0 ] ••• Rss[-N + 2 ] 


( 11 . 11 ) 


Rss[N- 1] Rss[N- 2] ••• R ss [ 0] 


Observe that R ss is Toeplitz, since all elements on each principal diagonal are equal. 
The correlation matrix is also Hermitian symmetric and positive semi-definite, but 
this is true regardless of the WSS assumption. We can form the correlation matrix 
from the reversed observation vector as: 


12 



R ss [0] JWl] "• R,s[N- 1] 

Rss[- 1] Rssl o] ••• Rss[N- 2] 

^ss [ — N + 1] i? ss [—iV + 2] ••• -R S s[0] 

Notice that because of the Hermitian symmetry property of the autocorrelation func¬ 
tion (i? ss [— l] — i?* s [Z]), it follows that 

R ss = K s (11.13) 

4. Cross-correlation Matrices 

For the purpose of representing the joint statistics of two random processes 
x[n] and y[n\, we define the cross-correlation function as: 

R X y[ni,n 0 ] = E{x[n 1 ] ■ y*[n 0 }} . (11.14) 

In the case of joint stationarity, the cross-correlation is also a function of the time lag 
l and the definition becomes 

.R I3/ [Z] = E {z[n] • y*[n — Z]} . (11.15) 

The cross-correlation vector between the signals s[n] and vector of observations x n is 
defined as: 

R SX [N - 1] 

r $x = E{s[n]-x* n }= \ (11.16) 

R sx [0] 

or in its time reversed form as: 

Rsx[ 0] 

f S i = E {s[n] • x*} = i . (11.17) 

Rsx[N-l] 
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More generally, given two sets of observations represented by the vectors x n of length 
N and y m of length M, we define their cross-correlation matrix as: 


R xy = £ {x n • y*J} = 


Rxy [0] 


R X y[ 1] 
R X y [0] 


Rzy[N-l] R xy [N- 2] 


R X y[—M + 1] 

R xy [—M + 2] 

Rxy [0] 


(11.18) 


The resulting matrix is not square (in general) and has no particular properties except 
for the repetition of elements along diagonals. 


B. DECIMATION AND UPSAMPLING 

The processes of decimation and upsampling are the basic tools of multirate 
digital signal processing. In the following sections we present two operators that 
extend sampling rate conversion as algebraic operations. 


1. Decimation Matrices 

The operations of decimation and upsampling can be represented by appro¬ 
priate matrix operations. Let 


S n 


s[n — 21V+1] s[n —21V + 2] 



(11.19) 


represent a sequence s[n], of length 2N ending with the observation s[n]. In that 
case, 


s[n — 2./V + 2] ... s[n + 2] s[n 


( 11 . 20 ) 


s[n — 2 N + 1] ... s[n + 3] s[n + 1] 
represent the even and odd sample observations as shown in Figure 2. The Decimation 
Matrix and Decimation Matrix with Time Delay are then defined as: 
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and 


D (0) = 


0 1 0 0 ••• 0 

0 0 0 1 ••• 0 


0 0 0 0 ••• 1 


D (1) = 


1 0 0 ••• 0 0 

0 0 1 ••• 0 0 


0 0 0 ••• 1 0 


and they are such that 


s e = D(°) 


S n 


and s Q = D^ 1 ) • s n . • 


( 11 . 21 ) 


( 11 . 22 ) 


(11.23) 


Both matrices D^ 0 ) and D- 1 ) are of size N x 21V. 


2 . 

Let 


Upsampling Matrices 


S n — 


s[n — N + 1] s[n — N + 2] 



(11.24) 


represent a sequence s[n] of length N ending with the observation s[n]. In this case 
the sequences 1 


0 s[n — N + 1] ... 0 s[n + 1] 0 s[n] 
s[n — IV+1] 0 ... s[n + l] 0 s[n] 0 


(11.25) 


1 Notice that superscripts are used in Equation 11.25, while subscripts are used in Equation 11.20. 
The reader should observe that s e ^ s e and s 0 ^ s°. 
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represent upsampling on s[n] up to the n th sample with no time delay or with time 
delay of 1 accordingly. The two upsampling operations can be represented by the 
Upsampling Matrix and Upsampling Matrix with Time Delay such that: 

s e = U<°> • s n , s° = UW • s n (11.26) 

with U< fc > G 7l 2NxN . The upsampling and the decimation matrices are related as 
follows: 

u (fc) = D (fc)r fc = 01 (H.27) 

Since upsampling followed by decimation has no effect on the original signal, we can 
easily show that 

D « . D (*)r = i N (H. 28 ) 

where I N is the identity matrix of size N. However, decimation followed by upsam¬ 
pling does not restore the original signal, so: 

U (fc) • U (fc > T ^ I 2N . (11.29) 

It is also easy to see that the reversal operation changes the unit delay introduced in 
D^, thus: 

£)(°) = D (i) and 5(i) = D (o) (11.30) 


C. CONVOLUTION MATRIX 

Consider the linear time-invariant operation where a signal x[n) is fed to an 
FIR filter with impulse response h[n] and of order P. The output y[n] will be the result 
of the convolution between x[n] and h[n] (y[n] = x [n] * h [n ]). Let us represent the 
sequence x[n) with an observation vector x of length N. We define the Convolution 
Matrix of h[n] associated with the vector x as the following (P + N — l)x N matrix: 
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The computation of an observation vector may end up significantly simpler 
with the use of decimation and upsampling matrices as well as convolution matrices. 
If we assume in particular the system displayed in Figure 3, we can represent the 
resulting observation vector x as a matrix product between the vector s and the 
corresponding matrices as: 

x = D ( °) • H - s. (11.32) 



Figure 3. Effect of Decimation and Filtering 

D. THE NOBLE IDENTITIES 

The Noble Identities [Ref. 6: pp. 119-120] represent one set of basic tools 
for polyphase representation for decimation and interpolation filters. They allow 
commuting of the decimation or the interpolation operator that leads to simplification 
of multirate structures. These identities are illustrated in Figure 4. In what follows, 
we try to form an equivalent of them with the use of vectors and matrices. Let us 
suppose that we have the finite impulse response of a Linear Time Invariant Filter of 
order P, h[n ] = {h[ 0], h[l ],..., h[P — 1]}. The ^-transform of this filter is denoted 
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H ( z ) and we already have the necessary tool for describing the convolution in terms 
of vectors and matrices from Equation 11.31. However, we need a way to describe 
H(z 2 ) in order to deal with the Noble Identities. The impulse response that would 
correspond to H(z 2 ) would be an “up sampled version” of h[n]. We can say that 
h^[n] = {[/i[0], 0, h[l], 0,..., h[P — 1], 0}, meaning that we have one non-zero sample 
every two samples. In this case, the corresponding convolution matrix is: 


H< 2 > = 


Mo] 

0 

0 

0 

0 

MO] 

0 

0 

Mi] 

0 

0 

0 

h[p ~ i] 

0 

0 

0 

0 

MP-1] 

0 

0 

0 

0 

MO] 

• 0 

0 

0 

0 

MO] 

0 

0 

••• h[P - 1] 

0 

0 

0 

0 

h[p - 1] 

0 

0 

0 

0 


N 


(11.33) 


The resulting matrix has dimensions (2 P 4- N — 1) x N. 

With the introduction of IT 2 ) we can obtain the the matrix form of the Noble 
Identities that correspond to Figure 4 as: 


y = H D (0) • x <=► y = D (0 > • H (2) • x 
as the first identity, and 


(11.34) 
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(11.35) 


y = D (0)T H x <*=*> y = H (2) • D (0)T • x 


as the second identity. 



Figure 4. The Noble Identities 
(a) Identity No.l (b) Identity No.2 

E. TIME VARYING SYSTEMS AND CYCLOSTATION¬ 
ARY PROCESSES 

In the classical theory of optimal filters we make the assumption that the 
process to be estimated is WSS and that it is sampled at the same rate as the 
observed process. As a consequence, optimal filters are also linear time-invariant 
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filters and the stationarity of the estimate is preserved. In the formulation of multirate 
systems however, this is not the case. Both decimation and upsampling are time- 
varying operations, and in some cases the output of a multirate system may exhibit 
periodicity. The simplest example is when a multirate system performs decimation 
and then upsampling by the same factor, say 2. In this case, if the input signal is 

x[n] = {x[0],x[l],x[2],...} (11.36) 

then the output signal is 

y[n) = {*[0],0, x[2], 0,x[4],0,...}. (11.37) 

The output y[n] is called a linear periodically time-varying (LPTV) signal [Ref. 6: 
pp. 130-132]. For the case where the input x[n] is random, wide-sense stationarity is 
lost in y[n ] due to the periodic appearance of the “zeros.” In this case, we need to 
introduce the concept of cyclostationary processes [Ref. 9]. A zero mean process x[n] 
is said to be wide-sense cyclostationary with period M if the correlation function is 
periodic such that 

R x [n + M, k + M] = R x [n, fc]. (11.38) 

In the example that was previously given, the cross-correlation between the input and 
the output of the multirate system will be: 

Rxy [^, 

It turns out that some of the correlation functions that are described in the follow¬ 
ing section pertain to cyclostationary processes. Furthermore, the period M of the 
processes is equal to the sampling rate factor. 

F. SECOND ORDER ANALYSIS OF LINEAR MULTI¬ 
RATE SYSTEMS 

In the development which follows, we address the effect of multirate systems 
on the correlation between different signals. Some typical problems are shown in 


= < 


R xx [n-k] , k = 0,2,4,.. 
0 , k = 1,3,5,.. 


(11.39) 
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Figure 5. In all cases under consideration, we wish to determine the cross-correlation 
matrix between the processes x{m] and y[n] given the fact the s[n] and y[n] are jointly 
stationary. The fact that we implement all the operations as linear algebraic products 
simplifies each problem considerably. The results are summarized below. 



B 
■1 

x[n] s[n] 

- > • - > 

y[n] 

a 

U 

x[tn] 

yW 

s[n] 

♦ - > 



1 

B 

S' w - 

x[m] s[n] 
- > + - > 

y[n] 


-J 

1 

x[m] 
-> 

y[n] 

s[n] 

« - > 

h[n] 

(0) 

- > 


s[n] 

t] 


(d) 

-j 


x[m] 

- y. 

y[n] 

(e) (f) 


Figure 5. Six Typical Cases of Cross-correlation 

(a) Filtering (b) Decimation 

(c) Filtering and Decimation (d) Decimation and Filtering 

(e) Filtering and Upsampling (f) Upsampling and Filtering 


1. Filtering ( Figure 5 (a) ) 

When the operation that is performed on s[n] is filtering with a FIR filter h[n] 
of order P, the resulting cross-correlation matrix can be computed as follows: 

R xy = E{x n - y* n T } = E {H s n • y*/} = H • R„. (11.40) 
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The matrix R iy has dimensions (P + N — 1) x jV, the same as H. 

2. Decimation ( Figure 5 (b) ) 

In the case of decimation by a factor of two, the resulting cross-correlation 
matrix can be computed as follows: 

RW = E {x m ■ y" + 4 = E {D<*> ■ s im+i ■ y£ +i } = D<*> • R,». (11.41) 

The corresponding cross-correlation function can be computed as follows: 

R^J[l] = E {x[m} ■ y*[2m — 21 + k]} — E {s[2m\ ■ y*[2m — 21 + &]} (11.42) 

which implies 

Rjg[t\ = R sy [2l - k}. (11.43) 

As it was previously stated, the cross-correlation function between x[m] and y[n], 
represents a cyclostationary process. The index k corresponds to the number of 
samples that occur within one period. Therefore, for a decimation factor of two, k is 
0 or 1. 

3. Combination of Filtering and Decimation ( Figure 5 
(c) & (d) ) 

The last two cases can be combined to use the effects of both filtering and 
decimation as in (c) and (d) in Figure 5. We have 

Rg> = E {x m ■ y£ +i } = E {D<*> • H ■ ■ y£ +t } (11.44) 

leading to 

r« = D «.H-R sj ,. (11.45) 

In the case of decimation prior to filtering, we can write 

RW = E {x m . yg, + 4 = E {H • D<‘> • s fc+1 ■ y£«} (11.46) 

and finally 

Rg> = H • • R ay . (11.47) 
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4. Combination of Filtering and Upsampling ( Figure 5 

(e) & (f) ) 

A similar case may occur when, instead of decimation we perform upsampling 
on the input signal. In the case where upsampling is performed prior to filtering, the 
resulting cross-correlation matrix is computed as: 

R<‘> = E {x„ • y’ r } = E {H ■ D<‘> r • s„ ■ y'/j (11.48) 

resulting in 

Rg> = H • D^ T • R, y . (11.49) 

When filtering is performed prior to upsampling, the resulting cross-correlation matrix 
becomes: 

R>? = E {x m • y-/} = E {D<‘> r H s.. y"} (11.50) 

which implies 

R(fc) = D (*)r . H . R y (11.51) 

A summary of the results from Equations 11.40 through 11.51 is provided in 
Table II. These results are fundamental to the derivation of an optimal filter that 
implements two sets of observations at different sampling rates. 
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Operation 

Result 

Filtering 

Riy — R" * R 5 J/ 

Decimation 

R<*> = D<‘> • R s , 

Filtering and Decimation 

Rg = D ( fc ) • H • R sy 

Decimation and Filtering 

Rg = H • DW • R sy 

Upsampling and Filtering 

Rg = H • . Rsy 

Filtering and Upsampling 

Rg = D( fc ) T • H • R sy 


Table II. Cross-correlation Table 
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III. THE MULTIRATE WIENER FILTER 


One major class of applications in multirate digital signal processing are the 
information fusion problems. The general concept in this class of problems is the 
following: Based on the combination of observations available from different sources, 
obtain the best possible estimate on a signal, to which the observations are somehow 
related. 

In accordance with the principle that was previously stated, we form the prob¬ 
lem that is illustrated in Figure 6. We wish to estimate a signal s[n], from which we 
can obtain two observations x[n] and y[m], related to the original signal, where the 
different indices denote different sampling rates. In particular, both signals s[n] and 
x[n\ are sampled at the highest possible rate, while y[m] is sampled at a lower rate; 
specifically the sampling rate for y[m] comes from decimation by an integer factor. 
We want to compute the estimate s[n] at the same rate as that of the original signal. 



Figure 6. Basic Multirate Optimal Filtering Problem 
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We consider the specific case where y[m] is at half the sampling rate of s[n] 
(decimation by a factor of two), and s[n] is zero-mean and wide sense stationary 
(WSS). However the problem can be easily extended to more than two observation 
sequences, each obtained with the use of a different decimation factor from the sam¬ 
pling rate of the original signal. 

We shall see that the estimate can be implemented as shown in Figure 7. The 
optimal filter consists of two branches that take as an input the observation signals 
x[n] and y[m]. After filtering from the corresponding functions and g^ k \ the high 
sampling rate output is decimated by a factor of two. Both outputs are added in 
order to provide the even and odd samples of the estimated signal s[n]. Finally both 
samples are clocked through a parallel-to-serial converter to produce the estimate. 
This structure is extendable to higher decimation factors. Thus, if the sampling rates 
of the observations were related to a decimation factor of three, the input signals to 
the parallel-to-serial stage would simply be three as well (s[3m+fc], where k = 0,1,2). 

A. THE MULTIRATE WIENER-HOPF EQUATION 

From the fact that we have two observation sequences available, we form two 
observation vectors which we call 'x.^m+k with N elements and y m with M elements. 
Since y[m] is sampled at half the rate, it is convenient to use the expression 2m + k, 
where k = 0,1 to represent the index of x[n]. The observation vector that is formed 
by combining x 2rn +fc and y m has the form: 
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Figure 7. Implementation of Multirate Wiener Filter 


x[2 m — N + 1 + k] 
x[2m — N + 2 + k] 

X-2m+k 

y m 

y[m — M + 2] 
y[m] 


x[2m + k] 
y[m — M + 1] 


(iiu) 


We start from the fact that the optimal mean-squared filter satisfies the orthogonality 
principle or, in other words, the error from the obtained estimate is orthogonal to the 
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observation vectors. Define e[n] as the error of the estimation, e[n] = s[n] — s[n], The 
orthogonality principle dictates that: 


t r i 

E < X2m+k • e*[n] ► = 0 (III.2) 

9m 

L J / 

where k = 0,1. 

We want to express the estimate s[n] as the result of linear filtering of the 
observation vector from a filter that obviously is not shift invariant. In addition, we 
can speculate that this filter is divided into two distinct filtering functions h[n,m ] 
(for a:[n]) and g[n,m] (for y[m}) and that they also depend on the index k due to 
decimation. We can write: 

N -1 M—l 

s[2m + k] = h^[2m + k,2m + k-l]- x[2m + k-l}+ g {k) [m,m-X]-y[m-X]. 

1=0 A =0 

(III.3) 

In order to make the notation more compact, let us express the impulse re¬ 
sponses of both filters as vectors. The estimate then becomes: 

s[2m + fc] = [ x£ m+fc |jd- h (III.4) 

L J g W 

where k = 0,1 or 

s[2m + k}= xl m+k ■ h w + 9m • g (fc) - (HI.5) 

Since reversing both terms in an inner product does not change the result, we can 
also use the following alternative expression: 

s[2m + k}= h (fc)T • x 2m+ fc + g (fc)T • y m . (HI.6) 

With the last expression in hand, we can now form the relations that satisfy the 
Orthogonality Principle. Equation III.2 becomes: 
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El 


^2 m+k 

9 m 


( s*[2m + k] — s*[2m + k ]) 


= 0 


and by substituting Equation III.6 we obtain 




^2 m+k 


( 

X2m+fc 

*T 

h( fc )* 

\\ 

E < 



• 

s*[2m + A:] — 


* 




{ 

9m 



9m 


g(fc)* 

)\ 


By making the necessary operations, 


~X-2m+k 




^2 m+k 


2^2 m+k 

*T 

* 

J 

i_ 


• s [2m + k] 

> = E < 

\ 


* 


• 


9m 

j 


{ 

9m 


9m 


g(fc)* 


we end up with 


•p* 

XXi xx 

* 

CsjT 

fff 


h( fc )* 


f * 

sx 

_ (R s ),r )* 

i 

of 


gW* 


~(fc)* 

sy 


where we have defined: 

~ E jxgm+fc ' 

R^ = £{ym-jC} 

R£>’ = E {i^ +k ■ y"} 


(III.7) 

(III-8) 

. (III.9) 

(111.10) 

(111.11) 

(III. 12) 

(III. 13) 
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r, x = E{s[2m + fc] -x^*} (III.14) 

r iy = E{s[8m + k\ -y* m }. (III.15) 

Note that R^) and are functions of k, while all other terms are not. Finally, by 
conjugating all terms in Equation III. 10, we obtain: 

r xx i Rg> i r h<*> 

Rif T I Ryy \ [ g (fc) 

for k = 0,1. 

We refer to Equation III. 16 as the multirate Wiener-Hopf equation. 

The filter that processes the estimate for s[n] is linear but time-varying. The 
vector that comes from the permutation of and gW has length M + TV, and 
since it is derived from the orthogonality principle, it produces the Minimum Mean 
Square Error (MMSE) estimate, given observations from the signals x[vi\ and y[m\. 
The estimate is given by Equation III.6 and implemented as the combination of two 
periodically time-varying linear filters, with a period equal to the decimation factor 
(2). To express the estimate in scalar form, we can write: 

N -1 M—1 

s[2m + k] = ^2 • z[2m + k - 1} + p (fc) [A] • y[m - X] (III. 17) 

1=0 A=0 

which is identical to Equation III.3, except we have observed that the solution to 
Equation III.16 does not depend on m, and have thus written the filter terms as 
h( fc )|7] and The optimal filter is then seen to have the implementation given 

in Figure 7. 
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We can use the fact that there is a solution to the filter to express the minimum 
error variance, as: 


< 7 g(k) = E {s[2m + fc] • £*[2m + k}} , k = 0,1. (HI-18) 

Note that the minimum mean square error is a function of the index k due to the 
periodically time-varying nature of the filter. By substituting the expression for 
e*[2m + k], we obtain: 


(Tg(k) = E < s[2m + k] • s[2m + k] — 


■*“2771-4-fc y m 


h( fc ) 


g 


(fc) 


\*1 


J J 


or: 


= E {s[2m + k ] • s*[2m + A;]} — E < 



\, , 

r 



h( fc )* 


E< 

\s[2m 4- k] • 

v *t 

A 2m4*fc 1 

1 V* T 

1 J m 

' 

g(fc)* 

► 

> 


= R ss [ 0] - 


fT I f(k)T 

sx I sy 


1 

,- 

* 

i_ 

- 

_ g (k) * _ 


R ss [ 0] - • hW* - f• gW 


(III-19) 


a 2 e (k) = R ss [ 0] - h.W* T • r„ - g (fc) * T - fg>. 


(III.20) 


B. COMPARISON OF SINGLE RATE AND MULTIRATE 
WIENER FILTERS 

We already have a form for the of, the MMSE that is the result of the multirate 
Wiener filter. It was previously stated that the derivation of the filter was based on 
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the same principles as the single rate Wiener filter. Therefore, it is only natural 
to pursue a quantitative comparison between the two. First we will compare the 
expression for the filter in these two cases and provide a connection between them. 
Secondly, if we denote as Oq the MMSE that is the result of the use of the single rate 
Wiener filter on only the high sampling rate observations x[n], we would like to show 
that (the MMSE from the multirate filter) is always a smaller quantity than <Tq 
and be able to derive an expression for the difference. 

To begin, let us consider the following: 


Lemma 1 Consider a square matrix K of the form 


K = 


A 

B 

b .t 

1 C . 


(III.21) 


Then the inverse K 1 can be expressed as 


K 1 = 


' A" 1 + G • E" 1 • G* T | 

i 

7 

a 

0 

i 

—E -1 • G* T 

E- 1 


where 


G = A -1 • B 
E = C — B* T ■ A -1 ■ B 


(111.22) 

(111.23) 


and where all of the cited inverse matrices are assumed to exist. 

Applying the results of Lemma 1 to the correlation matrix defined in 
left-hand side of Equation III. 16, we can write 

R -l _i P lrv —1 /~i*T I C' T? - 

xx + W ’ rhk ' 'Jfc I — Vjfc • -Cj/c 

— Efc _1 • Gf | Efc _1 

where 


(III.24) 



G„ = R„-‘ ■ R« 

E* =R js -RW".R„- 1 -R<‘ ) . 


(III.25) 
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Substituting in Equation III.24 to express the filter coefficients and g( k \ we 
obtain: 


h«0' 


✓ 

p -1 

1 o‘ 


—Gfc 


= < 




+ 


g<‘>_ 



0 

1 0 


I 


Et 


-Gf | I 


;(fc) 

• sy 

(II 


.26) 


where 0 denotes an all-zero matrix with appropriate dimensions and I denotes the 
identity matrix. In order to assure the existence of and g <k \ the matrix E fc must 
be invertible. We know that E*, is a submatrix that shares the main diagonal with 
R (fc) . In all cases of interest, where noise is present, is invertible. Since R^ is a 
positive definite matrix, E fc is also invertible and positive definite [Ref. 10: pp.333]. 
Therefore, we can express and g( fc ) as: 


h (fc) = R ia; _1 • i sx - G k ■ E*" 1 • (fg> - Gf - f SI ) (III.27) 

g (fe) = Er 1 • (fg> - Gf • f«) . (III.28) 

The term ho = Ru ' 1 • r SI is the well known solution of the single rate Wiener-Hopf 
Equation (see Equation II.2). Therefore, by combining Equations III.27 and III.28, 
we obtain: 


h (fc) = h 0 - G fc • g (fc) . (III.29) 

Equations III.27 through III.29 establish a relationship between the single rate Wiener 
filter and its multirate extension. This is particularly significant when we compute 
the mean squared error of the estimates to compare the two implementations, single 
rate and multirate. 
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The next step is to compute a new formula for of based on the previously 
stated equations. From Equation III.6, by substituting for the filter coefficients, we 
can write: 


s[2m + k] = x^ m+k ■ • g (fc > 

= ^*-(ho-G*.gW)+y'-g(*) 

(III.30) 

= xL+fc • ho - xL+fc ’ Gjk • g {k) + y T m ■ g (fe) 

= • ho + (Ym - * 2 m+k ■ Gfc) • g (fc) 

which implies 

s[2m + A:] = s 0 [2m + k] + (y T m - x£ m+fc • G fc ) • g (fc) . (III.31) 

where so[n] is the estimate from the Wiener filter based on the solution of the single 
rate Wiener-Hopf equations for x[n] as the only observation. We can see that the 
final estimate is the sum of so[n] plus a correction due to the new observations y[m). 

Now let us consider the mean-square error expression and show that the mul¬ 
tirate estimation gives a better estimate (in the mean squared sense) than the single 
rate Wiener filter with only one observation. By substituting Equation III.29 into 
Equation III.20, we obtain: 


a 2 e (k) =R SS [ 0] - h< fc > 


*T . £ 


ff (fc)*r . ~(/c) 

S *sy 


= - (h„ - G t • g (t ')' T f„ - g«* r • fW 

= S„(0] - h; 1 • f« . Gf ■ f« - gM" • fW 

- -«-' y 


(III.32) 
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or 


a li. k ) = (T o- S {k>T ■ (?S } - G* k T • r sx ) (III.33) 

where o\ = R ss [0] — ho T • r SI is the error covariance of the single rate Wiener filter. 
We can show that the second term in the right-hand side of Equation III.33 is always 
positive. Using Equation III.28, we obtain: 

g ( fc)*T . _ q*t . = (Er 1 (rj } - Gf • r sx ))* T • (fg) - Gf • fsi) 

= (*%> ~ Gf • f sx )* T • E," 1 • (fg) - Gf • r sx ) 

(him) 

where we have used the fact that Efc is Hermitian symmetric, so E£ T = Efc. Since Efc 
is positive definite, the last expression is a scalar quadratic form of a positive definite 
matrix ; therefore it is always positive. Since cr^(k) and Oq are both positive, we can 
conclude from Equation III.33 that 

< a\. (III.35) 

That is, the MSE of the multirate Wiener filter is always smaller than the MSE of 
the single rate Wiener filter. 

C. ESTIMATION OF FILTERED SIGNALS IN NOISE 

Based on the framework set in the previous section, we can address the follow¬ 
ing problem. A random signal s[n] which is zero-mean and WSS, is observed through 
two distinct channels at different sampling rates. We assume that a decimation by a 
factor of 2 takes place in one of the channels, while both observations are corrupted 
by noise of known covariance, independent on the signal s[n]. This is shown in Fig¬ 
ure 8, where the observations are denoted by x\ri\ and y[m], and the noise sequences 
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by u[n] and v[m) respectively. The goal is to compute an estimate for s[n], called 
s[n], that is the minimum mean-squared estimate for the desired signal, based on the 
two sequences of observations. The FIR filter that provides this estimate is called the 
multirate Wiener filter for the problem. 

In order to make the arguments easier to follow, we first consider a simplified 
case in which there are no transformation filters or, for both filtering functions in 
Figure 8, Q(z) = F(z) = 1. The last assumption makes the derivation of the filter 
less tedious. Following that original problem, we consider the generalization of the 
previous case with impulse responses 6{n) and 7 [n] of orders P and Q respectively. 



Figure 8. Block Diagram of Multirate Wiener Filter 

1. Unfiltered Observations 

We start with the simplified case, for which ©(z) = T(z) = 1 and assume that 
u[n) and v[m) as white noise sequences of variance and o* respectively. Based on 
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the problem parameters, we can compute each of the elements in both sides of the 
Wiener-Hopf equations for the multirate Wiener filter, as stated in Equation III.16. 
The resulting equations involve considerable simplifications, since both tt[n] and v[m] 
are white noise processes, uncorrelated with any other signal. The correlation matrix 
for the high sampling rate signal x[n] is: 

= E {x-2m+k ' = ^ + figm+fc) ‘ (s^m+A: ”h U 2m+k) } 

= E • S^ +i .| + E {sgm+fc • U^ +fe | . . 

, i , .1 (111.36) 

+E ju 2m+k ' S 2m+k } d” ^ y*-2m+k ‘ 

= r«") + ^.i n 

leading to 

R xx = R ( £M + al-I N (III.37) 

where the cross terms are zero and we have used the identity R ss = R* s (see Equa¬ 
tion 11.13), since the signal is stationary in the wide sense. The superscript (N, N) in 
the correlation matrix (R^ ,Ar )) represents the Toeplitz matrix produced from R ss [Z] 
and has N rows and N columns. However, the number of rows may not be necessarily 
the same with the number of columns, as it will be stated in the following derivations. 
The correlation matrix of the low sampling rate signal y[m] becomes: 

RJv = £{ym-y;'} = £{(D<‘’)-s 8m+i + v„)-(D( 0 )-s a , +l +v m )"'} 

= E {dm • g^ ■ sg, +t • DM'} + E {DM • s Sm+t • v"} 

+E {v„ ■ sZ+k ■ DM'} + E {v m • v"} 

= D (1J . r(2m,2M) . D (i)t + ^2 . j M 

(III.38) 
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which yields 

Ryy = D (1 ) • Rg M ’ 2M > • + a 2 v ■ 1 M (III.39) 

where E>(°) = DO) (see Equation 11.30). The product DO) • DO) r effectively divides 
the matrix dimensions by two. Actually, whether we use D^ 0 ) • D^ r or DO) • D0) r 
it is irrelevant for the derivation of Equation III.39, but we keep the first product 
for consistency. The fact that y[m] comes from decimating x[n) guarantees that it is 
WSS, since Ry y [l] = R xx [2l]. 

The cross-correlation term between the high and the low sampling rate signals 
is given by: 


RO)* = E {ism+k ■ y*J} = E {(s 2m+fc + u 2m+k ) • (D (fc ) • S 2m+k + v TO ) T ] 

= E {s 2m+ * • s£ +Jt • D<*> T } + E {s 2m+ Jt • 

+E {u 2m+Jfc • §Z +k • DO) T } + E (u 2m+fc • v*J} 


= RW2M) . £)(fc)r 


which implies 

RW = R<f- 2M) • D {k)T . 


(111.40) 

(111.41) 


We can also express the cross-correlation vectors, in the left-hand side of Equa¬ 
tion III. 16, in a similar fashion as: 


r sx = E {s[2m + k] ■ x * 2m+k } = E {s[2m + k] • ( s* 2m+k + u * 2m+k )} 

= E {s[2m + k] ■ s * 2m+k } + E {s[2m + k] ■ u 2m+k } (III.42) 

- ;W) 

A SS 
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and 


r sy = E {s[2m + k] • y^} = E { s[2m + k] • (D w • s* 2m+k + v m )} 

= D (fc ) • E { s[2m + fc] • s£ m+fc } + E {s[2m + k] • v^} (III.43) 


= fjw • 

Substituting Equations III.37 through III.43 into Equation III.16, we obtain: 



RW2M) . £j(fc)T 


h(fe) 


fWl) 

ss 

£><*> • 

| D (l) . R (2M,2M) . D (1)T + a 2 . 


g(fe) 


£)(*) • 


(III.44) 


Equation III.44 is the form of the Wiener-Hopf equation for the multirate Wiener 
filter when ©(z) = T(z) = 1. 

The structure of the correlation matrices R xx , Ry y and R!^ is as follows: 


R'rrrE — 


R ss { 0] R» s [- 1] 

R S s[ i] R S s[o] 

R SS [N- 1] R sa [N-2] 


while Ttyy is: 


Rss[—N + 1] 

R S s[~N + 2] 

Rssi 0] 


+<t,m 


■N 


(III.45) 




R ss [o] 

Rss{ 2] 


Rss[- 2] 

R ss [0] 


R ss [2M — 2] R SS [2M — 4] 


R SS [—2M + 2] 
R SS [—2M + 4] 

Rss[®} 


+ a 2 v -I M . (III-46) 
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Finally, has the form: 


i?ss[0 + A:] Rss {—2 + A:] • • * Rss [—2 (A/ — 1) + A:] 

Rss[l + k] R„[-l + k] ••• R„[-2(M-l) + l + k] 

R ss [N-l + k] R ss [N-2 + k] ••• R SS [N - 1 - 2(M - 1) + k) 

(III.47) 

Similarly, the structure of the cross-correlation vectors r SI and is as follows: 

R ss [ 0] 

Rssi 1] 

Rss{N - 1] 

and 

Rss[0 4- A:] 

7?ss[2 + A:] 

[2M - 2 + A:] 

2. Filtered Observations 

As previously mentioned, consider the original signal s[n] observed in both 
channels through FIR filters with impulse responses 9[n ] and 7 [n] with orders P 
and Q respectively. We also assume that the noise sequences u[n] and v[m] are not 
necessarily white. With the tools provided so far, we can obtain the correlation 
matrix of the combined signal observation vector. In the following, © and T denote 
the convolution matrices of 6[n ] and of 7 ( 71 ] respectively. Computing each term in 
Equation III. 16 separately, we obtain: 


(III.49) 




(III.48) 
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Rxx — E \*-2m+k • X-2m+k^ E {(©* • s 2m+k + U £m+k) ' (s*2m+k ' ® T + ^2m+fc)} 

= ®* • E {sg m+ fc • S|^ +A .| • © T + @* • E js 2m+k • ^■Sm+k'l 
+E |u 2m+k * Sg^ +fc | • © T + E {fig , n +k ’ ^ 2 m+fe} 

- ©* • R(f- N ') • © T + RM 

(III.50) 

which implies 

Rxx = © R^f ,JV,) • ©* T + Rf ’ N) (III.51) 

where the term N' = N — P + 1 denotes the dimensions that the matrix R ss should 
have in order for R zx to have dimensions N x N. The term ©* • § 2 m+ fc represents 
reverse filtering on the sequence s[n]. On the other hand, the term © • R ss ■ ©* T 
denotes both forward (zero-phase) and reverse filtering. 

Similarly the other diagonal block becomes: 

= «{y. ■ ff) = S {(6'" • 1“ ■ iw,+v.) • • r • D»'+<■;)} 

= e {dw • r- ■ § Sm+k ■ sK, +l ■ r r ■ dw} + e {d<'> • r* • i tm+k ■ <} 

+B {v m • sg, +t • r r • DW'} + E {v m • v"} 

= • r* • . r T • + r( m > m ) 

(III.52) 


which is written 


Ry y = D (1) • T • R(f'- w '). r* T • D (1)T + R(f M ) (III. 53 ) 
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where the term M' — 2 M — Q +1 denotes the dimensions that the matrix R ss should 
have in order for Rj^ to have dimensions M x M. 

The cross-term can be computed as: 


RW* = E {jtSm+k ■ Ym} 

= b {(©• • s Sm+t + ■ (b ( ‘< T ■ r" ■ s 2m+t + »;')' r } 

, , . (111.54) 

= ©* • E {sw • ■ r r • + ©• ■ e {*»,+* • v;u 

+E {iW* • s£ +1 } • r 1 . + E {u*„ +t . v"} 

= ©* • r^'’ m ') • r T • D( fc ) T 

which yields 

Rg = © • • T* T • D^ T . (III.55) 

Finally, 

R (fc)* T = £)(fc). r . ^ R *(v',A/')j * T . ©*T (III.56) 

By combining Equations III.51 through III.56, the left-hand side in Equation III. 16 
becomes: 


R (fc) - 


©•R(f>"'>-©* T + R™ 

f>w. r • (r *i N '’ M '^y T ■ ©* T 


© . . p*r . f )(k)T 


d^ 1 ) - r • R(f.^') • r* T • d (1 ^ t + r( m > m ) J 

(III. 57) 
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The cross-correlation vectors in the right-hand side of Equation III.16 become: 


r sx = E {s[n] • x* } = E {s[n] • (©* • s n + u„)*} 


— ®E {s[n] • s*} + E {s[n] • u* n } 


(III.58) 


which implies 


f„ = e fir.'). 


(III.59) 


Similarly 


fW = E {»[»] • yj,} = E {j[n] ■ ■ I" • s. + v„) } 




( 111 . 60 ) 


which becomes 


fg) = rn • r • - 1} . 


(III.61) 


We substitute Equations III.57, III.59 and III.61, in Equation III. 16 and obtain: 


R(*) • 


© • f(f- 1 ) 


f )( k ). r • 


(III.62) 


where is given by Equation III.57. The significance of Equation III.62 is that 
it provides a mean of computing the multirate Wiener filter coefficients and 
for k = 0,1, directly from the signal statistics (R ss ) and the filter impulse responses 
(© and T). 

A closer look at Equation III.62, reveals that part of it incorporates the single 
rate Wiener filter. Thus, we can produce a more general form for it, expressed from 


Equation III.63 that represents any signal that is subject to linear transformation (by 
6[n\) prior to the computation of the Wiener filter. 

(© R ss ■ ®* T + R„) • h = © • i ss (III.63) 

D. OTHER RELATED MULTIRATE PROBLEMS 

In the previous Chapter we proved that the estimate of the Multirate Wiener 
filter is optimal for the general problem illustrated in Figure 6. The next step will 
be to show how the multirate Wiener filter of Figure 8 can be modified and used 
in slightly different scenarios. A close examination of Figure 8 reveals two possible 
modifications. 

1. Perform decimation prior to filtering. The modified block diagram appears in 
Figure 9. 

2. Extend the result to decimation by a factor of K where K is an integer greater 
than or equal to 2. 

1. Decimation Prior to Filtering 

When the signal is decimated prior to filtering, also shown in Figure 9, the 
resulting observation vector y m is computed as: 

y m = r.D<°>s 2m+ ». (III. 64 ) 

The multirate Wiener-Hopf equation then has the explicit form: 

’ h« 1 [ ©-fif’ 1 ) 

R<*> • = (III.65) 

g (fc) r • fjw • f if' 


where 



R (fc) = 


© • R(f- vV,) • ©* T + R[ n ’ n) | © • R(f -' M>) • D^ T • r* T 

I 

r • D (fc) • R*( A/, . Ar '). @* t i d ( 1) • r • . r* T • d (1)t + r( a/ ’ a/ ) 

SS I SS V 

(III.66) 


and the dimensions are N' = N — P + 1 and M 1 = 2M — 2Q + 2. 



Figure 9. Block Diagram with Filtering Prior to Decimation 

2. Decimation by a Factor of K 

Let us return to the system of Figure 8 and suppose that the decimation is 
performed by a factor of K , where K >2. We can generalize the decimation matrix 
as D^, where the subscript K denotes the decimation factor and the superscript 
k = 0,1,2,..., K — 1 denotes the time delay that the matrix introduces. All matrices 
have dimensions M x K • M and can be computed as follows: 
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0 


D 


(fc) _ 
K ~ 


fe +1 


0,. 


K 


0 ••• 


fc +1 


0 


0 0 , --, 1 , 0,---,0 


K 


(III.67) 


All the identities that were valid for K = 2, can be generalized to any value K > 2 
as: 


U )c 


( fc ) = 


K 


(III.68) 


and 


fj£) = Df _1_A:) (III.69) 


where A; = 0, 1,2,... ,K — 1. The multirate Wiener-Hopf equation becomes: 


where 


R<*>- 


’ h« ’ 


e • f<r.» 

g(fc) 


r • £><*> ■ ?<"'■') 


(III.70) 


R (fc ) = 


© • Rif-*') • ©* T + Rif*) 


0 . R (n',m>) . d 


(k)T 

K 


f)W. p. p_*(Ar',jv'). 0 *t | D^f -1 ) • r • ,,M ') • r* T • + r[ m > m ) 

(III.71) 


and the dimensions are ./V 7 = AT — P + 1 and M' = 77 M — Q + 1. 
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The structure of the matrix equation depends on two factors. The correlation 
matrix has a 2 x 2 block structure that results from the presence of two observation 
sequences. In a more general situation of multiple observation sequences, the number 
of blocks would be equal to the number of those sequences. For instance, three 
observations would lead to a 3 x 3 structure of the correlation matrix and three 
distinct filtering functions, one for each observation sequence. On the other hand, 
the number of times that the filter has to be computed, depends on the number of 
samples between the observation on the high sampling rate and the observation on 
the low sampling rate. As an example, if we obtain three observations at decimation 
factors Di = 1 (no decimation), D 2 = 2 and D 3 = 4 of the original sampling rate, 
the multirate Wiener filter would involve three distinct filtering functions (say , 
and h^), and would have to be computed for k = 0,1,2,3 in order to cover the 
new observations from all signals. 


E. EXPERIMENTAL RESULTS 

The experiments were based on simulations created with the use of MATLAB. 
For the needs of the filter implementation, a function was created with the name 
MultiWiener.m . The source file is included in the Appendix, while the source code is 
included in an enclosed diskette. In all cases, the model that was used was a second 
order AR (auto-regressive) model with transfer function 1 /A(z) and generating noise 
variance cr^ = 1 as shown in Figure 10. The experiments were divided into two parts: 

1. The use of the function only in order to check its theoretical performance. In 
this case, the results included the filter coefficients and the computed minimum 
MSE, using different combinations of additive white noise. The performance 
was compared to that of the single rate Wiener filter. 

2. A set of two simulations where a signal was created, distorted and filtered ac¬ 
cording to the model. The first simulation was conducted without pre-filtering. 
For the second simulation, the pre-filtering transfer functions were obtained 
from the use of the Remez Algorithm. Both transfer functions described low- 
pass filters of different orders. The filter performance was compared to that 
of the single rate Wiener filter. 
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Figure 10. AR Model for Signal s[n 


1. Theoretical Performance 

For the first set of results without pre-filtering (@( 2 ) = 1 and r(.z) = 1), we 
considered several cases of different additive white noise variances <t^ and The 
filter orders for h (fc) and g (fc) remained the same (N = 6 and M = 4) for all cases. 
In addition, instead of presenting the filter vectors h w and g (k) , the following tables 
present their 2-norm (||h|| 2 ) as a metric for the weight of the filter coefficients. The 
results from this set are shown in Table III. 

A close look at Table III, provides some interesting observations. 

1. In the first case where both = 1 and = 1, h' 1 ' 1 has much larger weight 
than g (1) , since only the measurement from x[n\ is present. For k = 0 on 
the other hand, both filters have similar weights, since there are available 
measurements from both signals. 

2. In the second and third cases, where = 0 and o\ = 0 respectively, the 
weight of the filter coefficients is transfered completely to the respective filter 
(h^-* or g^), since the additive noise is zero. 
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Table III. Theoretical Multirate Filter Performance 

3. Finally in the fourth and fifth cases, where a 2 = 1000 and = 1000 respec¬ 
tively, the weight of the filter coefficients is also transfered by the most part 
according to what data has the smallest additive noise variance; that is why 
for both k = 0 and k = 1 the MSE is about the same when a 2 = 1000. The 
vahie of the second observation is evident in the last case, where the single 
rate Wiener filter has no effect and the resulting MSE is almost as much as 
the signal variance (/? ss [0]). 

2. Simulations 

The main purpose of the simulations was to confirm the theoretical predictions, 
by comparing the performance of both filters on simulated data. In both cases, the 
orders of the multirate Wiener filter were chosen as N = 12 and M = 8, while 
the order for the single rate Wiener filter was N + M = 12 + 8 = 20. Also, for 
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both simulations, the additive white noise variances were a\ = cr% = 1. In both 
simulations, the data generation was repeated 50 times. The following figures show 
the time averaged squared error for both the multirate and the single rate Wiener 
filter in each of the 50 experiments. 

The performance comparison for the first simulation, which was conducted 
without pre-filtering is shown in Figure 11. The theoretical minimum MSE from the 
multirate filter was a 2 e = 0.3959 for k = 0 and a] = 0.6116 for k - 1, yielding an 
average of 0.5038. On the other hand, the theoretical minimum MSE from the single 
rate Wiener filter was <j\ = 0.6572. The average computed from the 50 time averaged 
square errors shown in Figure 11, matches the <j\ values quite accurately. 

Multirate: N= 12, Single rate: N+M=20 
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Figure 11. Performance Comparison without Pre-filtering 

The second simulation was conducted with pre-filtering, using the Remez Al¬ 
gorithm in order to model the transfer functions. For the design of O(cj) the order 



52 



































was P = 6 and cut-off frequency has been chosen as 0.47T. For the design of r(u;) 
the order was Q = 12 and cut-off frequency 0.3 tt. Both selected frequency responses 
are presented in Figure 12. The performance comparison between the two filters is 
shown in Figure 13. The theoretical minimum MSE from the multirate filter was 
= 0.3801 for k = 0 and cv == 0.4093 for k = 1. The average of the two previ¬ 
ous values is = 0.3947. On the other hand, the theoretical minimum MSE from 
the single rate Wiener filter was cr| = 0.4665. The theoretical performance and the 
average from Figure 13 are very close. 



Figure 12. Magnitude Frequency Response for 0(u>) and F(u;) 

3. Some Notes on the Simulations 

The results from Figures 11 and 13 confirm the theoretical results, that is 
the Multirate Wiener filter performance is in all cases superior to that of the single 
rate Wiener filter. However, there are some points that should be pointed out about 
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Multirate: N= 12, Single rate: N+M= 20 



Figure 13. Performance Comparison with Pre-filtering 

the way the simulations were implemented. Operations like R„ = © • R ss • ©* T 
have the characteristic that they perform both ways, causal and anti-causal, allowing 
zero shift of the correlation function after filtering. On the other hand, operations 
like r sx = © • r ss cause a shift that turns out to be half the order of the filter (e.g. 
P/2). That shift was taken into account, both during the filter computation as well 
as during the performing of the simulation. The result was that both filtered versions 
of s[n], Sfl[n] and s 7 [n], were delayed by P/2 and Q/2 respectively. 
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IV. LINEAR PREDICTION FOR 
MULTIRATE SIGNALS 


In many applications, linear prediction is used as special type of optimal fil¬ 
tering. Although it is indeed an optimal filtering problem, linear prediction and the 
sister topic of autoregressive (AR) modeling are a very rich area with numerous spe¬ 
cific applications relating to signal processing. These include coding, spectral analysis, 
system modeling and others. Basically in linear prediction, the past values of a ran¬ 
dom process axe observed, and it is desired to estimate the current value (which is 
assumed unknown). Suppose that we form a vector out of the even and odd sample 
of a signal s[n] called a signal vector as: 

x[2n] 
x[2n — 1] 

The purpose of the analysis of this Chapter is to perform linear prediction on 
this kind of signal. Since this signal consists of two channels, the techniques used are 
multichannel techniques. We characterize a signal as multichannel, when possibly 
several correlated signals are received on different channels and then processed and 
analyzed together. A common example of a multichannel 2-D signal is a color image. 
The signal is two-dimensional but consists of three registered components representing 
red, green and blue intensities, or other combinations of values such as those used 
in various video standards. Although what was previously stated suggests otherwise, 
the motivation for applying multichannel linear prediction came from the Polyphase 
Representation. 

A. THE POLYPHASE REPRESENTATION 

The polyphase representation of a filter, is the foundation of a number of 
techniques in multirate signal processing. It permits great simplification of theo¬ 
retical results and also leads to computationally efficient implementations of dec- 



x[n] = 
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imation/interpolation filters, as well as filter banks, both single rate and multi¬ 
rate [Ref. 6]. To explain the basic idea, consider a filter with transfer function 
H(z ) = Y,n=-oo M n ] ’ z ~ n - By separating the even-numbered from the odd-numbered 
coefficients of h[n], we can write: 

H-oo +00 

H(z)= Y h[2n] ■ z~ 2n + z- 1 ■ Y h[2n + 1} ■ z~ 2n . (IV.2) 

n=—oo n=—oo 

Then by defining: 

-1-00 H-OO 

H 0 {z)= Y h[2n]-z- n , H x {z) = Y h[2n + l]-z~ n (IV.3) 

n=— oo n=~oo 

as the filters related to the even and odd samples of the impulse response respectively, 
we can write H{z) as: 


H{z) = H 0 (z 2 ) + z- 1 ■ ^ (z 2 ). (IV.4) 

The last representation holds whether H(z) is FIR or HR; causal or non-causal. It 
also provides an alternative representation of the filtering operation. Suppose that 
a deterministic signal x[n ] is filtered by H(z), and y[n\ is the output sequence. If 
we are interested in the decimated sequence y[2m] only, the polyphase representation 
becomes particularly efficient, as shown in Figure 14. Similarly, we can delay the 
input signal by one and use the same technique to obtain y[2m — 1] as shown in 
Figure 15. 

The Noble Identities, described in Chapter II, are at the basis of the results 
in Figures 14 and 15. The results that are shown in Figures 16 and 17 show that 
we can obtain the even or the odd samples of y[n] as a linear combination of filtered 
even and odd samples of r[n]. If we call X 0 (z) the ^-transform of x[2m] and X\ (z) 
the .Z-transform of a:[2m — 1], and Yo{z) and yi(z) the Z-transforms of y[2m] and 
y[2m — 1], we can form a vector from the two signals and represent the operations of 
both Figures 16 and 17 in matrix form as: 


56 




Figure 14. Even Samples y[2m] via the Polyphase Representation 



Figure 15. Odd Samples y[2m — 1] via the Polyphase Representation 
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Y„(z) 


Ho(z) 

1 ff.M 


X„(z) 

Yi M _ 


z~ l ■ H,(z) 

1 H„(z) _ 


x i(z) 


(IV.5) 


The matrix in Equation IV.5 is called the polyphase matrix. The importance of 
this decomposition comes from the fact that the transfer function H(z) is implemented 
at lower sampling rate. The structure of the polyphase matrix suggests two interesting 
problems: 


1. Devise a method to make the polyphase matrix diagonal, and decompose the 
even and odd components of y[n] into uncorrelated sequences. 

2. Obtain x Q [n] and Xi[n] from yo [n] and y\ [n]. In case the signal is stochastic, 
this decomposition would yield two innovations, from which the original signal 
can be reconstructed. 


B. MULTICHANNEL RANDOM PROCESSES 

A multichannel random process with M channels ( [Ref. 11], [Ref. 12], [Ref. 
13} ) is defined as a vector of M random processes: 


x[ 


n\ 


Xi[n) 

x 2 [n] 

x M [n] 


(IV.6) 


If all of these signals are jointly stationary, the multichannel autocorrelation function 
is a sequence ofMxM matrices, defined as: 


r x [Z] d = E {x[n] • x* r [n - /]} . (IV.7) 

For example, in the case M = 2, the ACF is given by 


58 




Figure 16. The Noble Identities for y[2m] 



Figure IT. The Noble Identities for y[2m — 1] 
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r x[^] — 


(IV.8) 


Rx\,X\ [^] Px\ y X2 [^] 


R-X2,X2 R] 


1. Multichannel Linear Prediction 

We can extend the concept of linear prediction to the multichannel case ( [Ref. 
14], [Ref. 13] ) and compute the estimate x[n] of a multichannel signal as a linear 
combination of P past values of the signal. The form of the estimate is: 


x[n] = -Af • x[n - 1] - A? • x[n - 2] - ... - Ap ■ x[n - P] (IV.9) 

where the coefficients are M x M matrices, defined as negated and Hermitian trans¬ 
poses for later convenience. The resulting prediction error e[n], is defined as: 

p 

e[n] = x[n] — x[n] = ^ A* T ■ x[n — *]. (IV.10) 

t=0 

For the estimate to be optimum in the mean square sense, the resulting mean squared 
error (MSE) defined as E {e[n) ■ e* T [n]} has to be minimized. The orthogonality 
principle in the multichannel case requires: 

E {x[n — /] • e* T [n]} = | / = ° (IV. 11) 

[ 0, l = 1,2,...,P. 

The matrix S e is called the prediction error covariance matrix and is not necessarily 
diagonal. 

2. The Multichannel AR Model 

Multichannel Unear prediction is closely related to multichannel auto-regressive 
(AR) modeling, where we assume that a multichannel signal x[n] can be represented 
by the following difference equation: 
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x[n] = —Aj T • x[n — 1] — A^ T • x[n — 2] — ... — Ap T • x[n — P] + w[n}. (IV.12) 


Here A* are the matrix coefficients and w[n] is a multichannel white noise process 
with correlation function 

r TO [Z] = E {w[n] • w* T [n - /]} = E w • S[l). (IV.13) 

Although the multichannel white noise process is uncorrelated with respect to the lag 
l, it may be correlated between channels. In other words, the covariance matrix E w 
is not necessarily diagonal. 

Linear prediction can also be performed in backward order. When we perform 
backward linear prediction, we compute the estimate x[n—P] of a multichannel signal 
as a linear combination of P future values of the signal. The form of the estimate is: 


x[n - P] = -Bf • x[n - P + 1] - Bf • x[n - P + 2] - ... - B> T • x[n]. (IV. 14) 

The AR model matrix coefficients or the linear prediction matrix coefficients can be 
found from the solution of the Normal Equations modified for the multichannel case. 
As in the single channel case, the Normal Equations can be obtained for forward or 
backward prediction. We will call the error covariance that results from the for¬ 
ward Normal Equations, and ££ the error covariance that results from the backward 
Normal Equations. If we call B, the coefficients for the backward prediction, we can 
write the two associated Normal Equations as: 
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(IV. 15) 


r x[0] 

r x[l] 

r x[-l] 

r x[0] 

rx[~P] 

r x [-P + l] 


for forward prediction, and 

r x[0] r x [-l] ... 

r x[l] r x [0] ••• r 

r x [P] r x [P — 1] 


r x[P] 


1 

»—1 

1 _ 


1 - 

M 

r x[ P - 1] 

• 

Ax 

= 

0 

r x[0] 


. Ap . 


-1 

• O 

_ i 


x[-P] 


I 


'sf 

-P+1] 

• 

Bx 

= 

0 

r x[0] 


Bp 


0 


(IV. 16) 


for backward prediction, where I denotes the identity matrix and 0 denotes the zero 
matrix, both of size M x M. From Equations IV. 15 and IV. 16 it can be seen that 
the error covariance matrices and E* are given by: 

p 

S e=E r *M-Ai (IV. 17) 

t=0 


p 

S e = E r *H]-B i (IV. 18) 

z=0 

where we define A 0 = Bo = I. 

3. Multichannel Levinson Algorithm 

The multichannel Levinson algorithm (Levinson-Wiggins-Robinson [Ref. 14], 
[Ref. 13] ) like its single channel counterpart, provides a recursive method for solving 
the multichannel Normal Equations. Instead of solving the normal equations directly, 
the algorithm computes the model parameters of order p from model parameters of 
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order p — 1. This is done for p = 0,1,2,..., P. Let us denote the model parameters 
of order p as: 



> 

_1 


>>' 

■> 

II 

1 

, Sp(— SO, Bp — 



A w 

- P 


i- 

Cd 

$ 

i_ 


, e;(= sj). (iv.i9) 


The algorithm is described by the following procedure: 


r x[l] r x [2] ••• r x [p] 


■ 


(IV.20) 


r p = (s;_ j- 1 • a : 




r; = (Sp_i) _I ■ Ap T 
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_i 
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Bp_i 


r -| 



0 
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— 

— 

• r p , Bp = 


— 

— 



Bp_i 





0 




0 


L J 


• r; (iv. 2 i) 


s p - Sp_i • (i r; • r p ), 

with initial conditions 


s^ = s' r . 1 -(i-r„-r p ) 


A„ = B„ = I , S„ = EJ = r^O] , Ai = r*[l]. (IV.22) 

The new matrices defined in the algorithm are the M x M forward and back¬ 
ward reflection coefficient matrices T p and T p and the partial correlation matrices A p 
and = A* T . Note that the matrices A p and B p have increasing dimensions. The 
algorithm involves block reversal , in the sense that the submatrices are reversed 
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in order within the matrix A p without being reversed internally. The algorithm does 
not present the symmetry of a single channel case. So in general, B p ^ A* T , Tp ^ T* r 
and Sp ^ E p . 

We can express the prediction error for the forward prediction as 

<*W = £Ai 0 ”" • x[n - i] (IV.23) 

i=0 

and for the backward prediction as 

4W = Y bI ? ) * T • x ( n - V + *]■ (IV.24) 

2=0 

With proper substitution in Equations IV.21, we can obtain a recursive form for the 
prediction error as: 

I -r; T 1 I" e p _ a [n] 

(IV.25) 

-T' p * T I 4-^ n ~ i] 

Equation IV.25 shows that the p th order forward and backward prediction errors can 
be computed from the corresponding (p — l) ih order errors. This is illustrated in 
Figure 18 where the p th order filter can be formed recursively from the (p — l) th order 
filter with the addition of a lattice 'section. 

This lattice representation is convenient since the parameters of earlier stages 
do not change as the order of the filter is increased. 

C. MULTIRATE LINEAR PREDICTION 

In this section we treat the decimated signal as a two-channel random process 
and we apply the multichannel Levinson algorithm. We start from obtaining a random 
process s[n] which is zero-mean and stationary in the wide sense. We can form 
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Figure 18. Lattice Section for Multichannel Linear Prediction 


a 2-channel signal by obtaining the even and odd samples separately. We define 
x 0 [n] = s[2n] as the process that comes from decimation by a factor of 2 on s[n] with 
zero time delay, and Xi[n] = s[2n — 1] as the process that comes from decimation by a 
factor of 2 on s[n] with time delay of one. We then define the following signal vector: 


x n = 


x 0 [n] 


s[2n] 

_ *i[n] _ 


s[2n — 1] 


(IV.26) 


The correlation function, can then be computed as: 
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r x [l] = E {x[n] • x* r [n — /]} 



s[2n] 
s[2n — 1] 


s*[2n — 2l] s*[2n — 1 — 21] 


= E{ 


s[2n] • s*[2n - 2/] 
s[2n - 1] • s* [2n - 21} 


s[2n] • s*[2n — 1 — 2/] 
s[2n - 1] • s*[2n -1-2 1] 


(IV.27) 


i? S5 [2Z] | R ss [2l + 1] 

Rss[2l-1] | Rss[2l} 


Observe that 


r _ -R„[2(-i)] I B„[2(-()+ l] 

B„[2(-() - 1] | fi„|2(-()] 

(IV.28) 

*;,[2 1 ] i b ;,[ 2 i - i ]1 , 

= = r" 1 . 

+ i «;,[2i] J 

Based on the definition of the multichannel Normal Equations (Equation IV. 15), 
the correlation matrix for forward prediction becomes: 
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Rss{ 0] 

Pss[l] 

••• P SS [2P] 

Rss[2P + l] 

Rss[~ i] 

Pss[0] 

••• Rss[2P — 1] 

Rss[2P] 

Rss[- 2] 

Rss[- 1] 

Rss[2P-2] 

Rss[2P- 1] 

jRss[-3] 

■Rss[ — 2] 

P ss [2P-3] 

Rss[2P-2] 

Rss[~2P] 

Rss[ — 2P +1] • • • 

Rss[o] 

Rss[l] 

:„[-2P-i] 

7M-2P] 

••• Rss{- 1] 

P 5 s[0] 


(IV.29) 


It is apparent that R/ is identical to the matrix that is used for single ch ann el linear 
prediction of order 2 P. On the other hand, the correlation matrix for backward 
prediction becomes: 


67 




-R ss [0] 

Rssi 1] 

••• i? w [-2P] 

R S si~2P + 1] 

Rss[- 1] 

■R S5 [0] 

••• R SS [—2P — 1] 

Rssi-2P] 

R ss [ 2] 


••• R $S [—2P + 2] 

i? ss [-2P + 3] 

Rss[ 1] 

Rssi 2] 

••• R SS [-2P + 1] 

R S s[—2P — 2] 

Rss[2P] 

R SS [2P + 1] ••• 

R ss [ 0] 

Rssi 1] 

-ss[2P - 1] 

Rss[2P] ,••• 

Rssi- 1] 

Rss[0} 


(IV.30) 

A closer look at Equations IV.26, IV.27, IV.29 and IV.30 reveals the fact 
that if we reverse the signal vector in Equation IV.26, the matrices R/ and R fe are 
interchanged. The last observation indicates that the choice for the order of the 
channels within the signal vector leads to equivalent results, a fact that is in line with 
what we should expect. 

1. Linear Prediction Output 

A possible way to check the output of linear prediction on the multirate pro¬ 
cess, is to try to form the prediction error. Based on Equation IV. 10, the prediction 
error in this case will be: 

p 

e[n] = x[n] — x[n] = ^ A* T ■ x[n — i] (IV.31) 

i=0 

or equvalently, 
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SoN 

£i[n] 


s[2n] 
s[2n — 1] 


+ £f=iAf- 


s[2n — 2i\ 
s[2n — 2z — 1] 


1 0 
0 1 


s[2n] 1 [ 4° 0) * 4 10) * 1 [ s[2n - 2] 1 

s[2n - 1] + 4 01) * 4 U) * s[2n — 3] 


(IV.32) 


s[2n] + 4°°)* • s[2n — 2} + 4 10 ' 1 * * s[2n — 3] + ... 
s[2n - 1] + 4 01) * • s[2n - 2] + 4 n) * • s[2n - 3] + ... 


Equation IV.32 provides an indication that at least one of the components of e[n] (in 
particular Si[n]) is essentially the same as the error e[n] that is the output of linear 
prediction in the single channel case on the same signal. This fact is also confirmed 
by the numerical results. 

2. Numerical Results 

For the purpose of performing multichannel linear prediction on the multirate 
random process, a MATLAB function was created with the name ChanLevin.m . The 
main purpose of this function was to implement the multirate Levinson recursion and 
show the results for A p , B p , r p , T' p , S p and E p for comparison and further analysis. 
Two autocorrelation functions were created for two AR models of orders two and 
three, (two and three poles respectively). The idea behind this was to exclude any 
chance that the observed experimental properties would be accidental. The first model 
was actually the same that was used in the Multirate filter experiments (shown in 
Figure 10); with generating noise variance = 1 and poles at 0.5 and 0.8. For the 
second model the generating noise variance was also cr^ — l and the poles were at 
0.4, 0.5 and 0.7. For both cases the linear prediction was performed for order P = 2. 
Since the underlying model is of order 2 and the theoretical autocorrelation function 
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was also computed for these results, linear prediction of higher order would not be 
necessary. 

The 2 nd order AR model is described by the difference equation 

s[n] = 1.3s[n — 1] — 0.4s[n — 2] + w[n] (IV.33) 

where a 2 = 1. This corresponds to a prediction error filter of the form A(z) = 
1 — 1.3z -1 T 0.4z~ 2 . The autocorrelation function terms for the process are: 


R ss [0] = 8.6420 R u [ 1] = 8.0247 R ss { 2] = 6.9753 
R ss [3} = 5.8580 R ss [4] = 4.8253 R ss [ 5} = 3.9297. 


(IV.34) 


Therefore, the correlation matrix blocks for the multichannel signal are given by: 


II 

o’ 

8.6420 

8.0247 

r x[l] = 

6.9753 

5.8580 


8.0247 

8.6420 


8.0247 

6.9753 

r x [2] = 

4.8253 

3.9297 





5.8580 

4.8253 





Finally, the results were: 


(IV.35) 




1.29 1.30 



-0.40 -0.52 



-0.52 -0.40 


r; 


1.30 1.29 

r 2 


0 0 

5 

. r ' 2 . 


0 0 



0 0 



0 0 


(IV.36) 
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1 

0 


1 

0 


0 

1 


0 

1 

— 

1.29 - 

-1.30 

, B 2 = 

0.40 0.52 


0.52 

0.40 


-1.30 -1.29 


0 

0 


0 

0 

- 

0 

0 


0 

0 


2.69 

1.30 

, S' 2 = 

1.00 1.30 



1.30 

1.00 


1.30 2.69 



(IV.37) 


(IV.38) 


The 3 rd order AR model is described by the difference equation: 

sjn] = 1.6s[n — 1] — 0.83s[n — 2] + 0.14s[n — 3] + w[ri\ (IV.39) 

where a 2 = 1. This corresponds to a prediction error filter of the form A(z) = 
1 — 1.6 z~ l + 0.83z -2 — 0.14z -3 . The autocorrelation fimction terms for the process 
axe: 


R ss { 0] = 13.1876 R ss [l] = 12.3347 R ss [ 2] = 10.5167 


(IV.40) 


R ss [ 3] = 8.4351 R ss [ 4] = 6.4942 R ss [ 5] = 4.8619. 


Therefore, the correlation matrix blocks for the multichannel signal are given by: 
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r x[0] = 


13.1876 12.3347 
12.3347 13.1876 


r x [2] = 


6.4942 

8.4351 


4.8619 

6.4942 


r x[l] = 


10.5167 

12.3347 


8.4351 

10.5167 


(IV.41) 


The obtained results were: 


ih 
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r; 
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, b 2 = 

0.8300 

1.1880 

1.1880 

0.8300 


-1.6000 

-1.7300 

— 
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— 

— 

-0.2240 

-0.1400 


0 
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0 


-0.1400 

-0.2240 


(IV.42) 


(IV.43) 



3.56 

1.60 


1.00 

1.60 

e 2 = 

1.60 

1.00 

, S' 2 = 

1.60 

3.56 


(IV.44) 


From Equations IV.36 through IV.44, we can make the following observations: 
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1. The error covariance for both forward and backward linear prediction is not 
diagonal. In addition, the diagonal elements are not equal. However, the 
cross-covariance elements are equal as they must be, since a covariance matrix 
is always symmetric. 

2. All 2x2 components are the reverses of their counterparts. In other words: 

Bf=A, w , r; = f p , S' p = t p . (IV.45) 

This makes sense since the signal vector is consisted of adjacent elements of 
the same random process. 

3. The observation made based on Equation IV.32 proved to be correct. The sec¬ 
ond column of A p (the one that corresponds to e\ [n]) has the same coefficients 
as the generating polynomial A(z). In addition, the first column of B p has 
also the same coefficients as the generating polynomial A(z) in reverse order. 

4. If we examine A p and B p columnwise, the number of non-zero elements are 
exactly the same as in the generating polynomial. Therefore, there is no 
computational loss or gain from the multichannel prediction. However, the 
number of matrix coefficients becomes half the selected order for large P. 

D. DIAGONALIZING THE ERROR COVARIANCE 

Despite the success in the implementation of multichannel linear prediction in 
the specific setup of even and odd samples of a random process and the remarkable 
symmetry that makes the particular use of multichannel prediction equivalent to its 
single channel counterpart, there is one problem that still remains unsolved. It is the 
fact that the error covariance still remains non-diagonal. If we recall the polyphase 
representation and the polyphase matrix from the beginning of the chapter, we would 
like to invert the model and design a whitening filter. In this way, we would obtain the 
innovation from the process by obtaining even and odd samples from two uncorrelated 
innovations. 

In this section we consider the application of suitable transformations that 
would diagonalize the error covariance matrix. We can consider two possible generic 
types of decompositions: 


73 



1. Eigenvalue decomposition: The eigenvalue problem is to determine the non¬ 
trivial solutions of the equation A • x = A • x. If A is an TV x TV matrix, the 
full solution of the problem includes N values that satisfy the equation called 
eigenvalues and arranged in a diagonal matrix A. The corresponding vectors of 
x are called right eigenvectors and can be arranged as the columns of a matrix 
V. If A is Hermitian symmetric (conjugate symmetric), the eigenvectors are 
always linearly independent and the eigenvector matrix V can be used to 
perform a similarity transformation, leading to the result V* r • A • V = A. 

2. Cholesky factorization: The Cholesky factorization [Ref. 10: pp.195] is another 
method that can be used to diagonalize a matrix. If A is Hermitian symmetric, 
it can be factorized as A = L-D-L* T , where L is lower triangular and invertible 
and D is diagonal. We can perform a similarity transformation on A, leading 
to the result L -1 • A • (L* T ) -1 = D. 

1. Diagonalizing with Eigenvalue Decomposition 

Since the error covariance matrix E^ is symmetric, it can be decomposed in 
terms as E{ = V • A • V* T , where V is the matrix of right eigenvectors and A is the 
diagonal matrix of eigenvalues. From Equation IV. 17 we can write: 

S' = V • A • V* T = £; r x [i] • A* (IV.46) 

i=0 

and therefore 

A = fv--r x [i].A t .V (IV.47) 

i=0 

where we have used the fact that V* T • V = V • V* T = I. If we define a new signal 
vector y[n] = V* T • x[n], then its correlation function becomes: 


T y [l] = E {y[n] • y* T [n — l]} = E {V* T • x[n] • x* T [n - l] ■ V} 


(IV.48) 


= V* T • r x [Z] • V. 


Furthermore, we can define the transformed coefficient matrices as 

Ci = V* T • Aj • V , * = 0,1,2,...,P 


(IV.49) 
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and write again Equation IV.47 as: 


A = Er y (i]-Ci. (IV.50) 

2=0 

This also means that the error signal vector for the prediction of y[n] can be written 
as: 


t' [n] = ES .0 Cf • y[n - i] = Ef,o V" • AT • V • V" ■ x[n - i] 
= V" • ELo A; t • x[n - i] = V* T • e[n). 


(IV.51) 


In conclusion, the transformation y[n] = V* T • x[n] yields the prediction ma¬ 
trices C i = V* T • A* • V, and a forward prediction error vector with uncorrelated 
channels. 

Let us now consider what happens with backward prediction. We have seen 
in the numerical examples (assuming that the signal on which linear prediction is 
performed comes from an AR model) that is the reverse of It can easily be 
shown therefore that the eigenvector matrices for and H b £ are V and V* T respec¬ 
tively. That means that we can determine a similar transformation for the backward 
prediction, defined as: y 6 [n] = V • x[n] with C b = V • A* • V* T . This transformation 
yields a backward prediction error with diagonal error covariance matrix. As a result 
of these considerations we see that we can apply a transformation to diagonalize ei¬ 
ther the forward or the backward error covariance but (unless V* T = V) we cannot 
diagonalize both simultaneously without changing the scale of the errors. 

2. Diagonalizing with Cholesky Factorization 

By an argument similar to the previous section, we can express using the 
Cholesky factorization as: 
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E' = L-D-L" = £>*[»]• A,. (IV.52) 

2=0 

By pre-multiplication and post-multiplication with L -1 and (L* T ) -1 respectively, we 
obtain 


D = ^ L" 1 • r x [i] • Aj • (L* T ) -1 . 
2=0 


(IV.53) 


We again define a new signal vector y[n] = L 1 • x[n], with correlation function: 


,[/] =E {y[n] • y* T [n - l}} = E {L^ • x[n] • x*> - /] • (L*T'} 


(IV.54) 


= L -1 • r x [/] • (L* T ) -1 

and by transforming the coefficient matrices to 

Q = L* T • Ai • (L* r ) -1 , i = 0,1,2,..., P (IV.55) 

we can write Equation IV.53 as: 


D = £>y[i]Cj. (IV.56) 

2=0 

Finally, the error signal vector for the prediction of y[n] can be written as: 


e'ln] = ZjL o Cf • y[n - i] = Zf= o L" 1 • A* T ■ L • L _1 • x[n - i] 

(IV.57) 

= • Ef=o A * r ■ x[n ~i}= L" 1 • e[nj. 

Thus, with the use of the transformation y[n] = L _1 • x[n], accompanied with 
Q = L* T ■ A i • (L* T ) _1 , we effectively make the signal vector produce a diagonal 
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error covariance matrix under forward linear prediction. Under backward predic¬ 
tion the case is similar to that of eigenvalue decomposition, and the corresponding 
transformation for the backward prediction is given by: y 6 [n] = (L* T ) -1 • x[n] with 
coefficient matrices C£ = L -1 • A, • L. The diagonal matrix D has all the diagonal 
elements equal. As in the previous case, the forward and backward predictions have 
two different diagonalizing transformations. 

3. Numerical Results 

In order to obtain results with the use of eigenvalue decomposition and Cholesky 
factorization, each of the 2x2 correlation function blocks (see Equation IV.41) were 
transformed according to the Equations IV.48 and IV.54. Afterwards, the multichan¬ 
nel linear prediction was performed with the function ChanLevin.m on the resulting 
correlation matrix. This time, only one AR model was used, the one that is shown in 
Figure 10 with generating noise variance = 1 and poles at 0.5 and 0.8. The linear 
prediction was performed for order P = 2. 

With the use of eigenvalue decomposition the results were: 

1.2325 0.4141 0.3115 0.0110 

1.4059 -0.3425 T[ 1.8310 0.5785 

— , — = — — (IV.58) 

o o r' 2 oo 

0 0 0 0 
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1 
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— 

— 

-1.2325 

-0.4141 

, B 2 = 

-0.3115 

-0.0110 

1.4059 

0.3425 


-1.8310 

-0.5785 

— 

— 


— 

— 

0 

0 


0 

0 

0 

0 


0 

0 



3.3955 

0 


2.4745 

1.4170 

e 2 = 

0 

0.2945 

, S' 2 = 

1.4170 

1.2155 


(IV.59) 


(IV.60) 


With the use of Cholesky factorization the results were as follows: 




1.0387 

1.6267 



0.2283 

-0.0185 



-0.1933 

-0.1487 


r; 


0.4833 

0.6617 

r 2 


0 

0 

> 

. r 2 . 


0 

0 



0 

0 



0 

0 


(IV.61) 
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1 
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1 

0 

0 

1 


0 

1 

-1.0387 

0.1933 

-1.6267 

0.1487 

, B 2 = 

-0.2283 

-0.4833 

0.0185 

-0.6617 

0 

0 


0 

0 

0 

0 


0 

0 

1.0000 

0.0000 

0.0000 

1.0000 

, S' 2 = 

0.3717 0.8167 

0.8167 4.4844 


(IV.62) 


(IV.63) 


From a close examination of Equations IV.58 through IV.63, we can make the 
following observations: 

1. Both transformations yield the expected results in terms of the error covariance 
E 2 for linear prediction. For the case of eigenvalue decomposition the error 
covariance matrix (Equation IV.60) was equal to the eigenvalue matrix. For 
the case of Cholesky factorization (Equation IV.63) the covariance of error was 
the identity matrix. 

2. There is no symmetry between the resulting 2x2 components in forward and 
backward prediction. As a consequence, E{ is different than E*, and E£ is 
not necessarily diagonal. 

3. Neither the elements of A p nor the elements of B p include the generating 
polynomial A(z), neither in the forward nor in the reverse order. 

4. Simulation Results 

A simulation was performed in order to check the performance of Unear predic¬ 
tion on data. The AR model that was used for data generation was the one shown in 
Figure 10 with generating noise variance <r% = 1 and poles at 0.5 and 0.8. The signal 
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s[n] formed two channels according to Equation IV.26 as x[n]. A linear predictor of 
order P = 2 was designed for x[n]. The residual error sequence e[n] was transformed 
by L -1 according to Equation IV.57 (Cholesky factorization), producing e'[n]. 

For both e^n] and [n] the autocorrelation fimction was computed, in order to 
check the theoretical prediction that they are orthogonal sequences. Also, their cross¬ 
correlation was computed in order to check the prediction that they are uncorrelated. 
The obtained results are presented in Figure 19, from we come to the conclusion that 
the implementation of Cholesky factorization transformed the two channels of the 
residual e'[n] to orthonormal sequences. 




Figure 19. Correlation Functions Between the Channels of e'fn] 

The residual e[n] was filtered from the inverse lattice based on the linear pre¬ 
diction of x[n] and expanded in order to check for perfect reconstruction of s[n]. The 
obtained estimate s[n] presented in Figure 20 is actually s[n] delayed by two samples, 
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since it comes from a two step ahead prediction. The figure confirms there is perfect 
signal reconstruction. 



Figure 20. Estimated Signal s[n] and s[n] from Inverse Lattice on e[ri 


The residual e'[n] was filtered through the inverse lattice, based on linear 
prediction on y[n] = L -1 • x[n] (Cholesky factorization), and expanded in order to 
check possible correlation with s[n]. The result, presented in Figure 21, shows that 
although the obtained estimate presents similar statistics, it is not an approximation 
of s[n}. 

We have shown that the vector e'[n] is stated as e'[n] = L -1 • e[n], therefore 
we can invert the situation and write e[n] = L • e'[n] . Since L is a lower triangular 
matrix, it can be written as: 

loo 0 

ho hi 


(IV.64) 
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Figure 21. Estimated Signal s[n] and s[n] Directly from Cholesky Factorization 


The significance of the Cholesky factorization approach comes from the fact that the 
elements of e[n] can be written as: 


eoM = loo ■ e'oH 

ei[n] = ho-£'o[n] + ln -e'lN 


(IV.65) 


with €q [n] and e'Jn] two orthogonal, uncorrelated sequences. As a consequence, given 
only e£)[n], we can reconstruct an approximation of the original signal from the inno¬ 
vation eoM- The additional information needed to reconstruct ei[n], and therefore 
the signal s[n], is the term e[ [n] which is imcorrelated to e' 0 [n]. This is described in 
Figure 22, which represents e' 0 [n ] and e' x [n] as two orthogonal (imcorrelated) vectors 
in a 2-D vector space. 
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Figure 22. Orthonormal Vector Representation of e'[n] 

The importance of the first channel of the residual e 0 [n] becomes evident from 
the fact that it can be used to produce an approximation of s[nj. We can use the 
principles of polyphase representation as they were stated in the first section of this 
Chapter and decompose the transfer function H(z ) of the AR model into Hq(z) and 
ifi(z), leading to: 


H(z) = 


(1 - 0.5 ■ z- 1 ) ■ (l - 0.8 ■ z- 1 ) 


(IV. 66) 


H 0 (z) = 


1 + (0.5 • 0.8) - z” 1 


(1 — 0.25 • z -1 ) • (1 — 0.64 • z~ l ) 


(IV. 67) 


H M (0-5 + 0-8) 

U ; (1 - 0.25 • z" 1 ) • (1 - 0.64 • z- 1 ) 


(IV.68) 


The first channel of the residual eo[n], was filtered from Hq(z) and the result was 
interpolated and compared with s[n]. The result is presented in Figure 23, where the 
estimate approximates the original signal s[n]. However, the same operation can be 
performed with e' Q [n) after multiplying it with the scalar quantity 1 0 q. 
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Figure 23. Signal s[n] and s[n] from Filtering and Interpolation on e 0 [n 


The importance of the last observation comes from the fact that we only need 
half of the samples of the residual to reconstruct an approximation of s[n], based on 
part of the prediction coefficients (the second column of A p as it appears in Equa¬ 
tion IV.37). In addition, the transformation of e[n] from L -1 produces two orthonor¬ 
mal sequences (e[)[n] and e'jfn]), from which we only need the first to obtain eo[n]. 
The last result provides to e' 0 [n] a functionality similar to that of the approximation 
coefficients in Wavelet decomposition (see [Ref. 1]). 
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V. CONCLUSIONS 


In this research we introduced multirate, multiresolution techniques to process 
random signals. In particular we focused on two objectives: estimation of a signal 
from a number of observation sequences at different sampling rates, and modeling of 
a signal presented in multiple channels formed from successive samples. The purpose 
of the latter was to design a set of filters to whiten the signed. 

For the estimation problem, the multirate Wiener filter was designed as an 
extension of the single rate Wiener filter with two observation sequences, one at 
high s am p ling rate and the other at low sampling rate. The technique is general 
enough to be extended to multiple observation sequences, each at a sampling rate 
that is a fraction of the basic sampling rate. The simulations that accompany the 
theoretical derivations, confirm the improved performance of the proposed technique. 
In addition, the proposed technique is extendable to fields such as Kalman filtering, 
adaptive filtering and also 2-D filtering and estimation. 

The signal modeling problem employed multichannel linear prediction. The 
motivation at the basis of the multichannel predictor is the representation of a signal 
by two systems and two uncorrelated white noise sequences. The two systems can be 
easily related either to the lower resolution and the higher resolution components of 
the signal, or to the even and odd samples of the signal itself. The fact that the two 
innovation sequences are white and uncorrelated is attractive in the sense that they 
represent the minimum amount of information (in the sense of minimum variance) 
required to reconstruct the signal. 

A possible application of this signal modeling is in speech coding using Linear 
Prediction, as shown in Figure 24. The overall information can be coded at a lower 
resolution, as an approximation of the original signal, and the details encoded in 
an uncorrelated innovation. The need to transmit the model parameters as side 
information represents a minimal overhead with respect to the overall amount of 
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data to be transmitted. 



Figure 24. Possible Application of Multirate Linear Prediction 

The main focus of this research has been on the development of the theoret¬ 
ical framework. Subsequent research should focus on the applicability of both the 
multirate Wiener filter as well as the multiresolution model to a diverse number of 
problems such as target tracking and speech coding. In addition this class of tech¬ 
niques is also attractive in image processing, as an extension to the multidimensional 
case. 
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APPENDIX. MATLAB CODES 


function [h,g,MSE]= MultiWienFIR(Rs,N,M,the,gam,varU,varV,i) 

% Function to find the time-variant FIR Wiener filter for a set of 
7% orders of N & M, given the ACF Rs[l] of a signal s[n]. The vectors 
7. ’the’ and ’gam’ denote FIR filters to which the signal is fed 
7 , before noise destortion. 

7. 

7. Syntax: [h,g,MSE] = MultiWienFIR(Rs,N,M,the,gam,varU,varV,i) ; where 
7. h : The part of the filter for x[n] (N x 1) 

7. g : The part of the filter for y[m] (M x 1) 

7 MSE : The minimum mean-square error 

7» Rs : The ACF from Rs[-N] to Rs[N] as a column vector. 

7. N,M : The orders for h,g respectively (N at least 1) 

7. the,gam : FIR filters impulse responses (P x 1) & (Q x 1) 

7. varU,varV : The additive noise variances of u[n] and v[m] 

7. i : Choise for even {0} or odd {1} time 

7. 

7. Note: If M= 0, we compute the classic Wiener filter 

7. An auxiliary function for Thesis. 

7. Dimitrios Koupatsiaris 

7. Last Revision: October 12, 2000 

k0= find(Rs == max(Rs) ) ; 

rsx= filter( the(end:-l:1),l,Rs ) ; ksx= find(rsx == max(rsx) ) ; 
Rxx= filtfilt(the,1,Rs) ; 

tempx= toeplitzC Rxx(k0:-l:k0-N+l),Rxx(k0:k0+N-l) . ’ ) ; 

RRxx= tempx + varU*eye(N) ; 

7. Reduce to classic FIR Wiener filter 
if M == 0, 

rr= rsx(ksx:ksx+N-l) ; 
h= inv(ERxx)*rr ; 
g= C] ; 

MSE= Rs(k0) - h’*rr ; 
return 

end 
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V. Decimation Matrix 


d= eye(M); D0= [] ; 

for k= 1:M, 

D0= [ DO , zeros(M,l) , d(:,k) ] ; 

end 

Dl= rev(DO) ; 

switch i 
case 0, 

D= DO ; 
case 1, 

D= rev(DO) ; 
otherwise, 

error(’Wrong choice for time parameter’) ; 

end 

rsy= filter( gam(end:-l:l) ,l,Rs ) ; ksy= find(rsy == max(rsy) ) 
Ryy= filtfilt(gam,l,Rs) ; 

tempy= toeplitz( Ryy(k0:-1 :kO-2*M+l) ,Ryy(kO:kO+2*M-l) . ’ ) ; 

RRyy = Dl*tempy*Dl’ + varV*eye(M) ; 

rxy= filter(the,1,Rs) ; Rxy= filter( gam,1,rxy(end:-l:1) ) ; 
Rxy= Rxy(end:-1:1) ; 

kxy= find( abs(Rxy) == max(abs(Rxy)) ) ; 

tempxy= toeplitz( Rxy(kxy:-1:kxy-N+l),Rxy(kxy:kxy+2*M-l).» ) ; 
RRxy= tempxy*rev(D’) ; 

RR= [ [RRxxjRRxy’],[RRxy;RRyy] ] ; 

rl= rsx(ksx:ksx+N-l) ; 

r2= rev(D)*rsy(ksy:ksy+2*M-l) ; 

rr= [ rl ; r2 ] ; 

H= inv(RR)*rr ; 

h= H(1:N) ; g= H(N+l:end) ; 

MSE= Rs(kO) - rl’*h - r2’*g ; 
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function [A,B,Sa,Sb,gam,gama]= ChanLevin(M,R) 

/, ChanLevin Computes the Levinson Recursion for an M-Channel signal 
/• given a ( M*P x M*P ) matrix R. The last matrix is a 

*/• block Toeplitz and represents the Autocorrelation Matrix 

7. for the M-channel signal. 

7» Syntax: [A,B,Sa,Sb,gam,gama]= ChanLevin (M,R) where 

7» A : The set of forward PEF coefficients (M*P x M) 

*/• B : The set of backward PEF coefficients (M*P x M) 

/• Sa : The forward covariance of prediction error (M x M) 

/• Sb : The backward covariance of prediction error (M x M) 

°/» gam : The reflection coefficients formed as a cell 

7. array of order P x 2 { GAMa,GAMb} 

7. gama : The reflection coefficients formed as matrix [GAMa.GAMb] 

7. Script file for THESIS work 
7. Advisors : 

7. Student : 

7. 

P= round(size(R,2)/M)-l ; 

A= eye (M); B= eye(M); 7. Initialize A & B 

Sa= R(1:M,1:M); Sb= Sa; 7. Initialize covariance of error 

gam= cell (P,2) ; gama= [] ; 


C.W. Therrien - R. Cristi 

Dimitrios Koupatsiaris 

Last Revision: November 03, 2000 


for i= 1:P, 

Amat= [] ; Bmat= [] ; temp= [] ; 
for j= 0:round(size(A,l)/M)-l, 

Amat= [Amat ; A(end-M*(j+1)+1 :end-M*j ,:)] 
Bmat= [Bmat ; B(end-M*(j+l)+l :end-M*j ,:)] 
temp= [ temp, R(1:M,M*(j+l)+l:M*(j+2)) ’ ] 

end 

Dp= temp*Amat ; 

GAMa= inv(Sb)*Dp ; 

GAMb= inv(Sa)*Dp > 
gam{i,l>= GAMa ; 


gam{i,2>= GAMb 


gama= [gama ; [GAMa,GAMb] ] ; 

A= [A;zeros(M)] - [zeros(M);Bmat]*GAMa 
B= [B;zeros(M)] - [zeros(M);Amat]*GAMb 
Sa= Sa*( eye(M) - GAMb*GAMa ) ; 

Sb= Sb*( eye(M) - GAMa*GAMb ) ; 


end 
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